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We show that seismic waiting time distributions in Cal- 
ifornia and Iceland have many features in common as, for 
example, a power-law decay with exponent a w 1.1 for inter- 
mediate and with exponent 7 ~ 0.6 for short waiting times. 
While the transition point between these two regimes scales 
proportionally with the size of the considered area, the full 
distribution is not universal and depends in a non-trivial 
way on the geological area under consideration and its size. 
This is due to the spatial distribution of epicenters which 
does not form a simple mono- fractal. Yet, the dependence 
of the waiting time distributions on the threshold magnitude 
seems to be universal. 



1. Introduction: Scaling laws in seismicity 

Earthquakes constitute an extremely complex spatio- 
temporal phenomenon, with the deformation and sudden 
rupture of some parts of the Earth's crust driven by con- 
vective motion in the mantle, and the radiation of energy 
in the form of seismic waves. Despite this complexity, cer- 
tain characteristics of seismicity can be captured by simple 
(empirical) self-similar laws. Thinking of complexity in the 
sense of the physics of complex systems (e.g. Gadomski 
et al. [2000]), one might argue that the observed scaling 
laws are in fact the hallmark of the crust being a com- 
plex system (see Mulargia and Geller [2003] for a recent 
review). For example, in a specified area covering many 
geological faults and over a given (large) time window the 
number of earthquakes N(m > mth) with a magnitude m 
larger than some threshold mth is given by the Gutenberg- 
Richter law (GR) [Gutenberg and Richter, 1949] which states 
that Iog 10 N(m > mth) oc ~b m t h with b ~ 1 [Frohlich and 
Davis, 1993]. Taking into account that the strain released 
during an earthquake is directly related to the moment M 
of the earthquake (by definition M = fj,A5 e where /1 is the 
shear modulus of the rock, A the area of the fault break, 
and 5 e is the mean displacement across the fault [Turcotte, 
1997]) which in turn increases exponentially with the mag- 
nitude (log 10 M — cm + d with c = 1.5, d = 10.73 [Stem 
and Wysession, 2002]), the probability distribution of the 
released strain turns out to be a power law, precisely the 
imprint of scale-free behavior. 

This self-similarity exhibited by GR, the independence 
of stress drop on magnitude [Ide and Beroza, 2001; Ka- 
gan, 2002], and recent laboratory measurements of coseis- 
mic slip resistance [Toro et al., 2004] are good evidence that 
the physics of earthquake rupture is the same for small and 
large earthquakes. It is still controversial if this is different 
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for the largest earthquakes (see [Kagan, 2002] for a discus- 
sion) because of the finite thickness of the brittle layer. 

Another example of a self-similar law in seismicity is 
Omori's law [Omori, 1894]. It states that immediately 
following a main shock there is a sequence of aftershocks 
whose frequency n decays with time t after the main shock 



n(t) = 



for t < tcutoff with p 



1 and fault- 



(K+t)P 

and magnitude-dependant constants K, C and t cuto ff [Utsu 
et al, 1995]. An obvious difficulty with Omori's law is the 
identification of aftershocks especially because aftershocks 
are not caused by a different relaxation mechanism than 
main shocks [Houghs and Jones, 1997; Helmstetter and Sor- 
nette, 2003]. Very recently, Baiesi and Paczuski [2004] pro- 
posed a solution based on a metric measuring the correlation 
between any two earthquakes. 

While the existence and limitations of GR and Omori's 
law are well tested and accepted, much less is known about 
other earthquake statistics and their (universal) properties. 
Here, we will focus on the distribution of time intervals be- 
tween successive earthquakes above a certain magnitude in 
a given area which is an important quantity characterizing 
earthquake occurrence. In the past, many different possi- 
bilities have been proposed for these waiting times, from 
totally random to periodic occurrence of large earthquakes. 
The most extended view is that of two separated processes, 
one for main shocks, which ought to follow a Poisson distri- 
bution [Gardner and Knopoff , 1974] (or not [Smalley et al., 
1987; Sornette and Knopoff, 1997; Wang and Kuo, 1998]), 
and an independent process to generate aftershocks. Very 
recently, Bak et al. [2002] suggested that waiting time dis- 
tributions can be described by a power-law with a cutoff for 
large waiting times which scales with the size of the consid- 
ered area and magnitude threshold, indicating a generalized 
scaling behavior. Corral [2003] further proposed that the 
"cutoff" is rather a transition between two different power 
laws. 

Here, we show that yet another power law regime exists 
for small waiting times. Our results also provide clear evi- 
dence that the additional transition point scales differently 
with area size. We attribute this to the generic structure 
of the epicenters' spatial distribution which is shown not 
to form a simple fractal but to have a rather complicated 
structure, possibly multifractal as proposed in [Hirabayashi 
et al, 1992; Goltz, 1998]. 



2. Waiting time distributions in California 
and Iceland 

To analyze the distribution of waiting times, we follow 
[Bak et al., 2002] and take the perspective of statistical 
physics by neglecting tectonic features and any classification 
of earthquakes as main shocks or aftershocks. We consider 
spatial areas and their subdivision into square cells of length 
L in kilometers. For each of these cells, only events with 
magnitude above a certain threshold mth are taken into ac- 
count. In this way, we obtain the waiting times n = ti+i — ti 
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in seconds between successive events at time ti and ti+i in 
each L 2 -cell. Combining the waiting times from all cells, the 
probability density function of the waiting times P mth ,L{T~) 
can be estimated. 

Analyzing data from southern California, Bak et al. [2002] 
proposed the following scaling law 

P mth ,L{T) = T- a j(rL d f/S (3 ) (1) 

with a ~ 1, d f ~ 1.2, f3 ~ 1, S = 10 mth and a scaling func- 
tion / depending only on the combined argument L df St. 
For small arguments, / is approximately constant not af- 
fecting the power-law (1/t) behavior; for large arguments, 
/ decays rapidly and in such a way that it seems to be 
consistent with a different power-law [Corral, 2003]. While 
there are good arguments that fi = b [Corral, 2003], a does 
not correspond to Omori's p as originally claimed by Bak 
et al. [2002] since the decay of the rate of aftershocks de- 
termines the corresponding all-return time distribution but 
not the waiting (or first-return) time distribution. It is an 
open question if df in Eq. 1 corresponds to the fractal di- 
mension D of the epicenter distribution: Corral [2003] has 
found Dq = 1.6 for southern California. It is also unclear 
if the analysis by Bak et al. [2002] has suffered from the in- 
completeness of their data for small magnitudes [see Wiemer 
and Wyss, 2000]. 

To test the general validity of the proposed scaling law 
(Eq. 1) and to clarify the scaling with L, we study 
two earthquake catalogues — one from southern Califor- 
nia (C, 19158 events) and one from southern Iceland (X, 
16286 events). For C, the coordinates of the polygon are 
(120.5° W,115.0°W) x (32.5°A, 36.0° N). Based on Wiemer 
and Wyss [2000], the reporting of events is assumed to be ho- 
mogeneous from January 1984 to December 2000 and com- 
plete at the level of m c = 2.4. For X, it was found that 
within (21.43° W, 19.8° W) x (63.62° A?", 64. 30° A'') the report- 
ing of events was homogeneous from July 1991 to December 
1995 and complete at the level of m c — 0.5 [Wyss, 2003]. 
The average waiting time differs by a factor of w 4 from X 
to C. 

For C, we find that indeed P mth ,L{T) decays as a power 
law with exponent a m 1.05 over some range of r. This can 
be deduced from Fig. 1 and Fig. 2 where Pm th ,L(T~) is plot- 
ted in terms of rescaled coordinates for different values of L. 
The z-axis is chosen as x = L d f S~ t with S = 10 mth , and 
the y-axis represents y — T a P mtht L(r). Thus, the constant 
regime in Fig. 1 and Fig. 2 corresponds to a power-law 
decay of P mth ,L(j) with exponent a. 

Yet, Fig. 1 and Fig. 2 also show that there is another, 
previously unnoticed power-law regime for smaller values 
of x: The power-law increase for the rescaled coordinates 
implies that the waiting time distributions decay with ex- 
ponent 7 « 0.64. Fig. 1 further shows that there is a rather 
sharp transition between the two power-law regimes char- 
acterized by a and 7, respectively. Moreover, for df — 1.0, 
all data collapse nicely onto a single well-defined curve for 
x < 10 6 and over 6 orders of magnitude including the tran- 
sition points. This data collapse implies in particular that 
the location of the transition point Ti between the a- and 7- 
regime depends on r, m t h and L only through the combined 
variable x with df = 1.0. For example, keeping m t h = 2.4 
fixed, we find Ti w 1.66h for L = 10km and Ti w 6.65h 
for L — 2.5km. Extrapolating to mth = 5, for instance, it 
follows Ti « 49.0h for L = 100km. 

However, the data do not collapse onto a single curve for 
x > 10 6 (see Fig. 1) implying that no single scaling func- 
tion / exists. Thus, the proposed scaling law given in Eq. 1 
cannot be strictly valid. This is further confirmed by vary- 
ing df. Already for df — 1.2 which was used in Bak et al. 
[2002] , the data do not collapse onto a single curve for small 



values of x. For df = Do = 1.6, the situation is even worse 
as Fig. 2 shows. Yet, with the exception of the curves for 
the smallest values of L, L = 2.5km and L — 5km, the data 
seem to collapse for x > 4x 10 and over 5 orders of magni- 
tude. Since the estimates of Pm th ,h(r) for large arguments 
and small values of L might suffer from insufficient length of 
the data set leading to an underestimation of the actual y 
values, a well-defined cutoff or transition point T2 between 
the power-law regime characterized by a and the regime for 
larger arguments showing a rapid decay (which might or 
might not be a power law) might exist. In contrast to Ti, 
T2 depends on the combined variable x with a different df, 
namely df = 1.6, and, thus, Ti and T2 scale differently with 
L. Consequently, differences in L are important and cannot 
be captured by a single exponent df. 

For X, we obtain similar results which are shown in Fig. 3 
together with data from C for comparison: For df = 1.0 
and fi = 0.95 and for x < 10 5 , the different curves collapse 
nicely over 5 orders of magnitude despite large variations 
in L and mth and their different geographical origin. In 
particular, within the statistical errors, there are no differ- 
ences between X and C in terms of a as well as 7. This 
suggests that a ~ 1.1 and 7 ~ 0.6 are universal. The data 
collapse further confirms that P mth ,L(j) depends on m t h 
only through the combined argument S^^r implying that a 
change in mth leads only to a rather simple rescaling of r. 
Besides, both for X and C we find f3 b. Thus, P mth ,L(r) 
scales with mth and according to GR: 

Pm th ,L(r) = P l (t /10» *"**). (2) 

Here, Pl does not only depend on L but also on the specific 
geographical area as Fig. 3 shows. Significant differences for 
large values of x between the curves with the same L but 
from different geographical regions are present. The waiting 
time distributions from C show a very sharp transition at 
and a much steeper decay for large arguments than the ones 
from X. These findings are evidence that even a modified 
version of Eq. 1 — incorporating the additional exponent 
7, for example — with a universal function / does not ex- 
ist. The data from X, similar to those from C, also do not 
collapse for large x despite identical mth- Thus, a simple 
scaling of Pm tK ,L(r) with L involving a single dimension df 
does generally not exist. As we show below, this is due to 
the fact that the distribution of epicenters does not form a 
simple fractal. 

3. Distribution of epicenters in California 

Any homogeneous (mono-)fractal is completely described 
by the capacity dimension Do alone. In the case of het- 
erogeneous (multi-)fractals, Do describes only one, albeit 
dominant, feature of the set as can be seen from the 
generalized definition of fractal dimension (the spectrum 
of generalized dimensions [Hentschel and Procaccia, 1983]) 

D q = lim^o (^-^fP) • Here, M q (L) = Y\ P«(L)« are 
the generalized qt\i moments of the probabilities Pi(L) = 
Ni(L)/N. Ni(L) is the number of events found in the ith 
cell of size L, and N is the total number. Do is recovered 
for q — 0. Two other prominent fractal dimensions, the in- 
formation and correlation dimension, result for q — » 1 and 
q = 2, respectively. Higher moments increasingly empha- 
size densely populated, i.e. seismically more active, areas 
of the set. D q may in principle be obtained for any value 
of q, especially desirable for negative q which would char- 
acterize the scaling properties of regions of low seismic ac- 
tivity. Yet, due to the relative sparsity of earthquake data 
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and its relative inaccuracy, we restrain our analysis to q > 
which will nevertheless allow us to judge the non-uniformity 
of the epicenter distribution. To obtain most reliable re- 
sults, we base our numerical approach on a generalization 
of the correlation-integral method as proposed by Pawelztk 
and Schuster [1987]. 

Figure 4 shows D q vs. q for C and for < q < 15. Inset 
is a double-logarithmic plot of M q (L) 1 ^ q ^ 1 ^ vs. L for the 
first three generalized dimensions and for q = 15, the high- 
est moment considered. D q is estimated from the respective 
slopes by fitting a straight line over the appropriate scal- 
ing region which ranges from about 1.5 to 180 km in this 
case, sufficiently large to believe in the fractal nature of the 
data. The scaling region holds well even for q = 15. We find 
Do = 1.60±0.13, in agreement with Corral [2003]. Further- 
more, Di = 1.38 ± 0.03, D 2 = 1.22 ± 0.05 and, effectively, 
Doo ~ 1. Doo is a measure for the fractal dimension of the 
most densely populated vicinities in the seismicity distribu- 
tion. Most notably, Do and Doo are identical within the 
statistical errors to the values of df for Tb and Ti, respec- 
tively. This suggests that P mth ,z,(r) is dominated by the 
most densely populated vicinities for small waiting times 
and the dominant feature of the spatial distribution of epi- 
centers captured by the fractal dimension Do only prevails 
for larger waiting times. This is exactly what is expected 
based on the associated difference in the local rates of seismic 
activity. Taking Doo — Do ~ 0.6 as a measure of inhomo- 
geneity, we conclude that the seismicity distribution in C is 
definitely heterogeneous and, thus, does not form a simple 
mono-fractal. For T, we reach the same conclusion. 

4. Discussion and conclusions 

The quality of our estimate of Pm th ,L(T~) depends on sev- 
eral aspects: first, on the completeness of the earthquake 
catalogue. Even if a catalogue is considered to be complete 
above m c , certain events are missing: Directly after a large 
earthquake many small events are lost in the seismic coda 
of the proceeding event. Hence, short waiting times will be 
underestimated and the error in n can be as large as the 
maximum "deadtime" after an earthquake. We, thus, have 
excluded waiting times less than one minute. 

Second, the magnitude of a given earthquake can vary as 
comparisons of different definitions of magnitude show. Yet, 
a detailed analysis of the data from I shows that P mth ,L{r) 
does not depend significantly on the definition of magnitude 
chosen. 

Third, while errors in the epicenter location can certainly 
lead to a wrong cell association and, thus, induce errors in 
the estimate of P mth ,L(r), the (implicit) assumption that 
earthquakes are a point process is much more severe. As 
shown by Kanamori and Anderson [1975], it is a good ap- 
proximation to relate the moment of an earthquake to the 
area of rupture by M a A z l 2 . Thus, A increases with the 
magnitude of the earthquake; for example, A 1 ^ 2 « 8km for 
m = 6 and A 1/2 w 80km for m = 8 [Turcotte, 1997]. Fortu- 
nately, for the vast majority of earthquakes A 1 ^ 2 is signifi- 
cantly less than the values of L studied here. This is espe- 
cially true for T since the rate of occurrence of the largest 
earthquakes is rather small, also preventing domination of 
the statistics by large events. 

We are, thus, confident that our estimate of Pm th ,L(r) 
is reliable and the scaling regime with exponent 7 for ar- 
guments smaller than Ti is not an artifact of our analysis. 
This is confirmed by the fact that the scaling of Ti and Tb 
with L can be directly related to the spectrum of general- 
ized dimensions, namely Doo and Do, respectively. While 
this directly implies that the proposed scaling law given in 
Eq. 1 is not valid, the comparison between vastly different 



scales and different geological areas suggests that a, 7 and 
Ti are universal. Our results further suggest that even if a 
power-law regime exists for arguments larger than T2 as pro- 
posed in [Corral, 2003], it might not be universal. It is up 
to future studies if the transition from 7 to a is associated 
with aftershock activity. 
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Figure 1. (color online) Rescaled waiting time distribu- 
tions for California with a = 1.05, df = 1.0, S = 10 m '\ 
rrith = 2.4, j3 = 0.95. The solid line corresponds to a 
power-law decay of P mth ,L{r) with exponent 7 = 0.64. 




Figure 2. (color online) Rescaled waiting time distribu- 
tion for California as in Fig. 1 but with df — 1.6. 




Figure 3. (color online) Combined rescaled waiting time 
distributions for Iceland (solid symbols, mth = 0.5) and 
California (open symbols, mth = 2.4) with a = 1.075, 
P = 0.95, d f = 1.6, S = 10 mth . Since the average error 
in epicenter locations is smaller for Iceland (< 1km) than 
for California (f» 1.7km), we have used different values 
for L. 




Figure 4. Spectrum of generalized dimensions D q for 
the epicenter distribution of California. Inset: general- 
ized moments, from which D q is obtained. 



